Filling a Silo with a Mixture of Grains: 
Friction-Induced Segregation 
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We study the filling process of a two-dimensional silo with inelastic particles by simulation of 
a granular media lattice gas (GMLG) model. We calculate the surface shape and flow profiles 
for a monodisperse system and we introduce a novel generalization of the GMLG model for a 
binary mixture of particles of different friction properties where, for the first time, we measure the 
segregation process on the surface. The results are in good agreement with a recent theory, and we 
explain the observed small deviations by the nonuniform velocity profile. 
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Mixtures of grains tend to separate as a response to vir- 
tually any type of external perturbation Q| , an effect that 
can be a major problem in some practical situations, and 
useful in others. Indeed, understanding the underlying 
processes of segregation phenomena in granular mixtures 
is an intriguing problem of interest to scientists from a 
wide range of disciplines. 

Segregation occurs when a mixture is poured onto a 
horizontal base and a pile builds up. When the grains of 
the mixture differ in size, the large grains tend to gather 
at the bottom of the pile , while when the grains have 
different shapes, the more faceted grains are found pref- 
erentially near the top. An even more surprising effect 
can be observed when a mixture of small rounded and 
large faceted grains is poured between two parallel plates: 
Stratification is observed — i.e., the grains organize spon- 
taneously into stripes ||| . 

Particle size distribution crucially influences segrega- 
tion. However, the particles generally differ not only in 
size but also in other properties, such as frictional ones. 
Frictional effects play a relevant role in band formation 
in 3-D rotating drums, in segregation in thin rotating 
drums, and also segregation and stratification in 2-D si- 
los. Thus the study of the particular case where segrega- 
tion is caused solely by differences in friction coefficients 
tells the contribution of friction to these segregation pro- 
cesses. 

The description of segregation in the case of a bidis- 
perse pile is nontrivial. A theory was proposed for the 
two-dimensional case E9j, which is a generalization of a 
method developed H for avalanches in piles of monodis- 
perse particles that treats the static bulk and the flu- 
idized surface (the rolling phase) separately, and a set of 



continuum equations describes the dynamics of the flow- 
ing region and the interactions between the two phases. 
Solutions have been found for the steady state filling of 
a 2D silo for the case of monodisperse particles. The 
theory provides also predictions for the surface and seg- 
regation profiles in a binary mixture consisting of grains 
with different friction properties but the same size |Q . 

It is important to test these theoretical predictions 
for two reasons. First, the complete description requires 
a number of constitutive relations between the relevant 
variables of the problem; these are usually unknown and 
are assumed to have analytically treatable functional 
forms, which should be verified. Second, in order to de- 
rive closed formulae for the profiles, several assumptions 
are made, so the limits of these approximations are of in- 
terest. Hence we have carried out a program of computer 
simulations using the granular media lattice gas (GMLG) 
model, which has been used for several granular systems, 
such as pipe flows, shaken boxes, and static piles l^)- 
An advantage of the simulation approach is that quanti- 
ties can be measured that are inaccessible in laboratory 
experiments. 

In the GMLG model §, one generalizes a fully discrete 
hydrodynamic algorithm Q in order to include energy 
dissipation through particle collisions and friction. The 
indistinguishable point-like particles are either at rest or 
else they travel with unit momentum along the bonds of a 
triangular lattice. The particles are scattered at the lat- 
tice nodes at integer times and then they are transferred 
to the nearest-neighbor sites in parallel. 

We adapt the GMLG algorithm to the case of two types 
of particles with different friction properties, which we 
call up and down particles 0). Using probability vari- 
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ables, we introduce material parameters in a stochastic 
way. The restitution coefficient is described by a param- 
eter e, which is the probability that energy is conserved 
in a collision (an example of the application of this rule 
is shown in Fig. |l]a). Momentum is conserved if the par- 
ticles in a collision are not connected, even in an indirect 
way, to the wall. 

The compact static part of the pile behaves like a 
solid with a large mass, where friction effects are taken 
into account. When moving particles interact with the 
bulk, their momentum can be transferred through the 
force chains to the walls of the vessel. For two types of 
grains, we define four different friction coefficients tx a p, 
(a, j3 Gt, I), giving the probability of a moving particle of 
type a to stop when arriving at a bulk site containing a 
particle of type (3 (Fig. |l|b). Bulk particles are rest parti- 
cles that are supported by another bulk particle or by the 
vessel. The particles have equal size, their distinction be- 
ing introduced through different friction coefficients 0] . 

We study the steady state filling process of a two- 
dimensional "silo", i.e., a long rectangular box of lateral 
size L. The silo is filled with a steady flux of particles Q 
next to the right wall. Two typical snapshots are shown 
in Fig. H for a monodisperse system of particles and for 
a binary mixture. The theory focuses on the limit of 
very slow, but still continuous filling. In general, the 
steady state slope depends on the incoming particle flux 
but this dependence vanishes at low rates. We find that 
for Q < 0.5 p d a a '^° ep , the slope will indeed be indepen- 
dent of the filling rate, while we still observe smooth and 
non-intermittent growth. For even smaller incoming flux 
(Q < 0.1 in the above units) the pile grows intermit- 
tently. This avalanche regime was studied in Ref. || in 
detail. 

In steady filling of a monodisperse species, theory pre- 
dicts that the thickness of the rolling particle layer, R(x), 
decreases linearly down the slope j4j, and that the local 
slope 9(x) is close to the angle of repose 9 rep everywhere 
but near the bottom of the pile (see also Ref. M for the 
case of a sandpile in an open cell), 

R(x) = %x, 6{x) = 9 rep - —. (1) 

vL jx 

Here 7 describes the rate of exchange between the bulk 
and the rolling phase, and v/j is of the order of the grain 
size. We assume the flow velocity v to be constant in 
space and time. 

Next we test these predictions by calculating R(x), de- 
fined in the simulation by the average number of rolling 
particles at each x, and the height of the pile, dh(x)/dx = 
tan #(2;). We take data when the steady state is reached, 
and we repeat the measurements at time intervals dur- 
ing which the pile grows about two lattice sites in height. 
We average the data over an ensemble of 10 — 100 sys- 
tems. The number of particles at the end of each run is 
typically of the order of 10 5 . 



Figure |3J shows the simulated profiles for the monodis- 
perse case. Figure |^a illustrates the fact that [see Eq. (|])] 

R(x) oc x (2) 

with reasonable precision. The slight deviation from lin- 
earity can be understood due to corrections to the ve- 
locity profile. Equation (^|) comes from a conservation 
of grains argument, assuming that the velocity of the 
rolling grains v is constant along the slope. However, at 
the bottom, where the slope is less steep we expect the 
particles to slow down and this is also observed in the 
simulations. If we take this into account theoretically 
by introducing a first-order correction to the velocity, 
v(6) = vq + X(9 — 9 rep ), we arrive at an implicit for- 
mula for R(x) (solutions for A = were derived in Ref. 
I): 

R t (x) ee R(x) + In (r(x) V A=-^-x. (3) 

By transforming our data with this solution, we obtain 
an excellent fit (Fig. ||a inset) . 

From the height profile (Fig. ||), we see that the slope 
at the upper part is almost constant, and a region can 
be found at the bottom where the surface flattens out. 
The singularity at x = [Eq. (Q)] cannot be seen, since 
both in real systems and in simulations there is a cutoff 
due to the grain size, but the profile does bend up (for 
small e) near x = 0. This effect has also been observed 
experimentally ||ic|| . The singularity can be avoided by 
re-defining the interaction term between the rolling and 
static grains. 

If the silo is filled with a mixture of particles (Fig. |]b), 
the growth process is considerably more complicated. In 
general, instead of one single angle of repose there are 
two continuous sets of generalized angles of repose for 
the two species || - denoted by 9] and #j - since the lo- 
cal critical angles depend also on the volume fraction of 
the species in the bulk. The curves are characterized by 
four variables 9 a p , which are the critical angles for parti- 
cles of type a rolling on a pure static phase consisting of 
grain type 0. As in Ref. [Q, we assume that both curves 
are constant. Then the number of critical angles reduces 
to two constants 9f = #n = #T1 an< ^ @l = @U = %t ( or 
equivalently, /Ltf and fj,±) |3j. 

The rolling phase in the steady state is described by 
the number of rolling particles for each species R a (x) 
(with a G|, I), and the static phase by the bulk volume 
fractions $f (x), and <&i(x), where $j + $j = 1. We 
distinguish two regions, where different analytic results 
have been obtained Q . The outer region includes almost 
the entire pile surface except for a narrow zone, the in- 
ner region, close to the bottom of the pile. We will focus 
on the flow properties in the outer region. The profile 
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R(x) = R-\(x) + Ri(x) is given by Eq. (1) while the pro- 
files of the components can be expressed in terms of an 
exponent r: 



i? T (x) : 



R(x) 



i+St(#)' 

R(x) 



(4) 
(5) 



The exponent r plays also a role in the determination of 
the bulk volume fractions $| and 



$ T (x) = 



#l(z) 



(6) 
(7) 



Here we assume a homogeneous rolling phase with con- 
stant velocity v = v-\ = and Qf and Qi are the fluxes 
of each species. The exponent r depends on the structure 
of the collision matrix describing the interaction between 
the bulk and rolling phase. 

To analyze the simulation data, we calculate the expo- 
nent r from the measured R(x) and R a (x) profiles for all 
x. The most significant test of the theory is the existence 
of r; if r is well-defined, the volume fraction profiles can 
be calculated and compared to the measured ones. 

Figure ||b shows a snapshot of the simulation for a 
mixture with ^ = 0.25 (0 T = 57°) and fj, l = 0.152 
(6 1 = 48°). We see the segregation of the mixture with 
the more sticky species (f) at the top of the pile. Figure 
|]a shows typical moving particle profiles. It is apparent 
that there is a slight higher order deviation from linear- 
ity in the case of R(x), as opposed to the prediction of 
Eq. (|^). The discrepancy is small, but it is significant 
enough that the theoretical profiles based on a linear ap- 
proximation of the curve do not fit the simulation results. 
However, if we use the measured R(x) for calculating the 
rolling particle profiles, Eqs. (^) and (^) hold to a better 
approximation (this bias will be discussed below). 

Most crucial is to verify Eqs. (|4|) and (||) by calcu- 
lating r. We present results for two sets of critical an- 
gles A and B, with ipA ~ 0.15 and ipB ~ 0.5, where 
ip = 0f — 6 1. Figure || demonstrates that the exponent is 
well defined in both cases except for a region at the top 
of the pile. The measured exponents are ta — 0.19 ±0.02 
and vb = 0.51 ± 0.03. This result is reassuring since r 
is expected to be of the order of ip Q. The exponent 
slightly depends on the Qi/Q-\ ratio, but is independent 
of the total flux provided Q is sufficiently small. 

With the help of the exponent we can obtain the vol- 
ume fraction profiles based on Eqs. (||) and (Q), and 
compare them to the measured ones (Fig. |]b). We find 
good agreement, except for a slight deviation in the inner 
region. We find stronger segregation for larger ip (or r); 



for the small angle difference ipA the segregation of the 
species is less pronounced. 

Although the numerical results fit the continuum the- 
ory, we observe some deviations. At the top of the pile, 
we see a discrepancy both at the rolling particle profiles 
and when calculating the r exponent. Here the dynam- 
ics is significantly different from what is considered in 
the continuum model: moving particles tend to be in 
free flight after collisions with the pile surface. 

Relation (|J) is not satisfied rigorously, as mentioned 
above. The reason is that, similar to the monodis- 
perse case, the x component of the velocity is not uni- 
form. For a weighted sum of the rolling particle densi- 
ties, however, linearity should hold to a better approxi- 
mation: v^(x)R^(x) + vi(x)Ri(x) oc x. We show a jus- 
tification of this ansatz in Fig. |4ja, by using the mea- 
sured v^(x) and v^(x) functions. The functional forms 
of these profiles can be approximated by vq — —, where 



0.9 ± 0.01 
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update step 



and a = 1.6 ± 0.2 



lattice unit 
update step 



w update step update step 

are fitting parameters. The functional form of the veloc- 
ities is consistent with the correction we found for the 
monodisperse system v a (x) = vq + X(9 — 9 a ), a s|, |, as- 
suming that the angle of the pile behaves approximately 
as Eq. (0). 

In general, we expect deviations from the theoretical 
predictions as an increasing number of particles are in 
free flight above the pile. If many particles tend to get 
detached from the surface due to elastic collisions such a 
contribution should also be incorporated into the theory. 
In the numerical results presented above, both particle- 
particle and particle-wall collisions are almost perfectly 
inelastic (the coefficient of restitution is around 0.25). 
Test runs for a more elastic medium show that the expo- 
nent describing the R a (x) profiles is no longer constant 
as a function of x indicating the limits of the theory. 

In summary, we have simulated the filling process of 
a two-dimensional silo, and our results for inelastic par- 
ticles compare well with the predictions of recent theo- 
ries of surface flow of granular mixtures. We find small 
corrections that are accounted for by a modified set of 
equations. Thus, we have verified the main assumptions 
of the theory, pointed out the limits of the approxima- 
tions involved, suggested improvements, and found good 
agreement between the simulation results and the im- 
proved theory. Future work is needed to test further the 
limitations of the theory and to generalize the model to 
more complex situations like polydispersity in the parti- 
cle size distribution. 

This work was supported by OTKA (T016568, 
T024004), MAKA (93b-352), the Heinrich-Hertz 
Stiftung, and the NSF. 
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FIG. Q (a) Illustration of the implementation of the restitution coefficient through a triple collision. If energy is 
conserved — with probability e — the particles are scattered. If the collision is dissipative — with probability (1 — e) — the 
particles stop, (b) An example of the application of the friction rule, when a moving particle of type a (white) arrives 
at a bulk site with a particle of type (3 (shadow). With probability /J, a p, the particle of type a loses its energy and 
becomes part of the bulk, while with probability 1 — ji a ^ its momentum is conserved and the node is not considered 
to belong any longer to the static phase. 

FIG.|| Two simulation snapshots for the steady filling of a silo, which is being filled with (a) uniform particles (b) 
a mixture of two different types of particles. Note that each pixel represents a lattice node that can be occupied by 
up to six particles. In (b) black and white dots denote lattice nodes where the J, and f particles are in majority, 
respectively (0f > 6^). The more sticky j particles are found preferentially at the top. 

FIG. | Profiles of the growing pile in the steady state of the monodisperse system, (a) The average number of rolling 
particles, R{x). The inset shows that Rt(x) [see Eq. (||)] is indeed proportional to x to a good approximation, with 
fitting parameter A = 0-1 ^date step ' 3 ) The ne ight of the surface, h(x). At each measurement the mean height is 
subtracted, since the pile is constantly growing. On Figs|||^ quantities are given in natural units of the simulation 
method: Lengths are measured in lattice constants, times in update steps. 

FIG. | Profiles in case of a granular mixture for tps ~ 0.5. (a) Rolling particles. The continuous lines are fits 
calculated using Eqs. (Q) and (||). The inset shows that v^R-^ + V[R[ oc x. (b) Volume fraction of the [ particles in 
the static bulk. (For the f particles = 1 — $j.) The continuous line shows the calculated profile according to Eq. 
(0) using the measured r exponent. 

FIG.|| The value of the exponent r calculated at each site x for ipA = 0.15 and ipB — 0.5; r is well-defined except 
for the uppermost part of the pile. 
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